"\n"
"// Returns net particle density to the atmosphere top\n"
"//\n"
"//                    | Zenith   .\'\n"
"//                    | angle  .\'\n"
"//                    |      .\'\n"
"//                    |    .\'\n"
"//                    |  .\'\n"
"//           _ _ _ _ _|.\'\n"
"//            A       |\n"
"//  fAltitude |       |\n"
"//           _V_ _ _ _|\n"
"//              .  \'  \'  \'  .\n"
"//           .\'               \'.\n"
"//          \'                   \'\n"
"float2 GetNetParticleDensity(in float fAltitude,            // Altitude (height above the sea level) of the starting point\n"
"                             in float fCosZenithAngle,      // Cosine of the zenith angle\n"
"                             in float fAtmBottomAltitude,   // Altitude of the bottom atmosphere boundary\n"
"                             in float fAtmAltitudeRangeInv  // 1 / (fAtmTopAltitude - fAtmBottomAltitude)\n"
"                             )\n"
"{\n"
"    float fNormalizedAltitude = (fAltitude - fAtmBottomAltitude) * fAtmAltitudeRangeInv;\n"
"    return g_tex2DOccludedNetDensityToAtmTop.SampleLevel(g_tex2DOccludedNetDensityToAtmTop_sampler, float2(fNormalizedAltitude, fCosZenithAngle*0.5+0.5), 0);\n"
"}\n"
"\n"
"void ApplyPhaseFunctions(inout float3 f3RayleighInscattering,\n"
"                         inout float3 f3MieInscattering,\n"
"                         in float cosTheta)\n"
"{\n"
"    f3RayleighInscattering *= g_MediaParams.f4AngularRayleighSctrCoeff.rgb * (1.0 + cosTheta*cosTheta);\n"
"    \n"
"    // Apply Cornette-Shanks phase function (see Nishita et al. 93):\n"
"    // F(theta) = 1/(4*PI) * 3*(1-g^2) / (2*(2+g^2)) * (1+cos^2(theta)) / (1 + g^2 - 2g*cos(theta))^(3/2)\n"
"    // f4CS_g = ( 3*(1-g^2) / (2*(2+g^2)), 1+g^2, -2g, 1 )\n"
"    float fDenom = rsqrt( dot(g_MediaParams.f4CS_g.yz, float2(1.0, cosTheta)) ); // 1 / (1 + g^2 - 2g*cos(theta))^(1/2)\n"
"    float fCornettePhaseFunc = g_MediaParams.f4CS_g.x * (fDenom*fDenom*fDenom) * (1.0 + cosTheta*cosTheta);\n"
"    f3MieInscattering *= g_MediaParams.f4AngularMieSctrCoeff.rgb * fCornettePhaseFunc;\n"
"}\n"
"\n"
"// This function computes atmospheric properties in the given point\n"
"void GetAtmosphereProperties(in float3  f3Pos,\n"
"                             in float3  f3EarthCentre,\n"
"                             in float   fEarthRadius,\n"
"                             in float   fAtmBottomAltitude,\n"
"                             in float   fAtmAltitudeRangeInv,\n"
"							 in float4  f4ParticleScaleHeight,\n"
"                             in float3  f3DirOnLight,\n"
"                             out float2 f2ParticleDensity,\n"
"                             out float2 f2NetParticleDensityToAtmTop)\n"
"{\n"
"    // Calculate the point height above the SPHERICAL Earth surface:\n"
"    float3 f3EarthCentreToPointDir = f3Pos - f3EarthCentre;\n"
"    float fDistToEarthCentre = length(f3EarthCentreToPointDir);\n"
"    f3EarthCentreToPointDir /= fDistToEarthCentre;\n"
"    float fHeightAboveSurface = fDistToEarthCentre - fEarthRadius;\n"
"\n"
"    f2ParticleDensity = exp( -fHeightAboveSurface * f4ParticleScaleHeight.zw );\n"
"\n"
"    // Get net particle density from the integration point to the top of the atmosphere:\n"
"    float fCosSunZenithAngleForCurrPoint = dot( f3EarthCentreToPointDir, f3DirOnLight );\n"
"    f2NetParticleDensityToAtmTop = GetNetParticleDensity(fHeightAboveSurface, fCosSunZenithAngleForCurrPoint, fAtmBottomAltitude, fAtmAltitudeRangeInv);\n"
"}\n"
"\n"
"// This function computes differential inscattering for the given particle densities \n"
"// (without applying phase functions)\n"
"void ComputePointDiffInsctr(in float2 f2ParticleDensityInCurrPoint,\n"
"                            in float2 f2NetParticleDensityFromCam,\n"
"                            in float2 f2NetParticleDensityToAtmTop,\n"
"                            out float3 f3DRlghInsctr,\n"
"                            out float3 f3DMieInsctr)\n"
"{\n"
"    // Compute total particle density from the top of the atmosphere through the integraion point to camera\n"
"    float2 f2TotalParticleDensity = f2NetParticleDensityFromCam + f2NetParticleDensityToAtmTop;\n"
"        \n"
"    // Get optical depth\n"
"    float3 f3TotalRlghOpticalDepth = g_MediaParams.f4RayleighExtinctionCoeff.rgb * f2TotalParticleDensity.x;\n"
"    float3 f3TotalMieOpticalDepth  = g_MediaParams.f4MieExtinctionCoeff.rgb      * f2TotalParticleDensity.y;\n"
"        \n"
"    // And total extinction for the current integration point:\n"
"    float3 f3TotalExtinction = exp( -(f3TotalRlghOpticalDepth + f3TotalMieOpticalDepth) );\n"
"\n"
"    f3DRlghInsctr = f2ParticleDensityInCurrPoint.x * f3TotalExtinction;\n"
"    f3DMieInsctr  = f2ParticleDensityInCurrPoint.y * f3TotalExtinction; \n"
"}\n"
"\n"
"void ComputeInsctrIntegral(in float3    f3RayStart,\n"
"                           in float3    f3RayEnd,\n"
"                           in float3    f3EarthCentre,\n"
"                           in float     fEarthRadius,\n"
"                           in float     fAtmBottomAltitude,\n"
"                           in float     fAtmAltitudeRangeInv,\n"
"						   in float4    f4ParticleScaleHeight,\n"
"                           in float3    f3DirOnLight,\n"
"                           in uint      uiNumSteps,\n"
"                           inout float2 f2NetParticleDensityFromCam,\n"
"                           inout float3 f3RayleighInscattering,\n"
"                           inout float3 f3MieInscattering)\n"
"{\n"
"    // Evaluate the integrand at the starting point\n"
"    float2 f2PrevParticleDensity = float2(0.0, 0.0);\n"
"    float2 f2NetParticleDensityToAtmTop;\n"
"    GetAtmosphereProperties(f3RayStart,\n"
"                            f3EarthCentre,\n"
"                            fEarthRadius,\n"
"                            fAtmBottomAltitude,\n"
"                            fAtmAltitudeRangeInv,\n"
"                            f4ParticleScaleHeight,\n"
"                            f3DirOnLight,\n"
"                            f2PrevParticleDensity,\n"
"                            f2NetParticleDensityToAtmTop);\n"
"\n"
"    float3 f3PrevDiffRInsctr = float3(0.0, 0.0, 0.0);\n"
"    float3 f3PrevDiffMInsctr = float3(0.0, 0.0, 0.0);\n"
"    ComputePointDiffInsctr(f2PrevParticleDensity, f2NetParticleDensityFromCam, f2NetParticleDensityToAtmTop, f3PrevDiffRInsctr, f3PrevDiffMInsctr);\n"
"\n"
"    float fRayLen = length(f3RayEnd - f3RayStart);\n"
"\n"
"    // We want to place more samples when the starting point is close to the surface,\n"
"    // but for high altitudes linear distribution works better.\n"
"    float fStartAltitude = length(f3RayStart - f3EarthCentre) - fEarthRadius;\n"
"    float pwr = lerp(2.0, 1.0, saturate((fStartAltitude - fAtmBottomAltitude) * fAtmAltitudeRangeInv));\n"
"\n"
"    float fPrevSampleDist = 0.0;\n"
"    for (uint uiSampleNum = 1u; uiSampleNum <= uiNumSteps; ++uiSampleNum)\n"
"    {\n"
"        // Evaluate the function at the end of each section and compute the area of a trapezoid\n"
"        \n"
"        // We need to place more samples closer to the start point and fewer samples farther away.\n"
"        // I tried to use more scientific approach (see the bottom of the file), however\n"
"        // it did not work out because uniform media assumption is inapplicable.\n"
"        // So instead we will use ad-hoc approach (power function) that works quite well.\n"
"        float r = pow(float(uiSampleNum) / float(uiNumSteps), pwr);\n"
"        float3 f3CurrPos = lerp(f3RayStart, f3RayEnd, r);\n"
"        float fCurrSampleDist = fRayLen * r;\n"
"        float fStepLen = fCurrSampleDist - fPrevSampleDist;\n"
"        fPrevSampleDist = fCurrSampleDist;\n"
"\n"
"        float2 f2ParticleDensity;\n"
"        GetAtmosphereProperties(f3CurrPos,\n"
"                                f3EarthCentre,\n"
"                                fEarthRadius,\n"
"                                fAtmBottomAltitude,\n"
"                                fAtmAltitudeRangeInv,\n"
"                                f4ParticleScaleHeight,\n"
"                                f3DirOnLight,\n"
"                                f2ParticleDensity,\n"
"                                f2NetParticleDensityToAtmTop);\n"
"\n"
"        // Accumulate net particle density from the camera to the integration point:\n"
"        f2NetParticleDensityFromCam += (f2PrevParticleDensity + f2ParticleDensity) * (fStepLen / 2.0);\n"
"        f2PrevParticleDensity = f2ParticleDensity;\n"
"\n"
"        float3 f3DRlghInsctr, f3DMieInsctr;\n"
"        ComputePointDiffInsctr(f2ParticleDensity, f2NetParticleDensityFromCam, f2NetParticleDensityToAtmTop, f3DRlghInsctr, f3DMieInsctr);\n"
"\n"
"        f3RayleighInscattering += (f3DRlghInsctr + f3PrevDiffRInsctr) * (fStepLen / 2.0);\n"
"        f3MieInscattering      += (f3DMieInsctr  + f3PrevDiffMInsctr) * (fStepLen / 2.0);\n"
"\n"
"        f3PrevDiffRInsctr = f3DRlghInsctr;\n"
"        f3PrevDiffMInsctr = f3DMieInsctr;\n"
"    }\n"
"}\n"
"\n"
"\n"
"void IntegrateUnshadowedInscattering(in float3   f3RayStart, \n"
"                                     in float3   f3RayEnd,\n"
"                                     in float3   f3ViewDir,\n"
"                                     in float3   f3EarthCentre,\n"
"                                     in float    fEarthRadius,\n"
"                                     in float    fAtmBottomAltitude,\n"
"                                     in float    fAtmAltitudeRangeInv,\n"
"							         in float4   f4ParticleScaleHeight,\n"
"                                     in float3   f3DirOnLight,\n"
"                                     in uint     uiNumSteps,\n"
"                                     out float3  f3Inscattering,\n"
"                                     out float3  f3Extinction)\n"
"{\n"
"    float2 f2NetParticleDensityFromCam = float2(0.0, 0.0);\n"
"    float3 f3RayleighInscattering = float3(0.0, 0.0, 0.0);\n"
"    float3 f3MieInscattering = float3(0.0, 0.0, 0.0);\n"
"    ComputeInsctrIntegral( f3RayStart,\n"
"                           f3RayEnd,\n"
"                           f3EarthCentre,\n"
"                           fEarthRadius,\n"
"                           fAtmBottomAltitude,\n"
"                           fAtmAltitudeRangeInv,\n"
"                           f4ParticleScaleHeight,\n"
"                           f3DirOnLight,\n"
"                           uiNumSteps,\n"
"                           f2NetParticleDensityFromCam,\n"
"                           f3RayleighInscattering,\n"
"                           f3MieInscattering);\n"
"\n"
"    float3 f3TotalRlghOpticalDepth = g_MediaParams.f4RayleighExtinctionCoeff.rgb * f2NetParticleDensityFromCam.x;\n"
"    float3 f3TotalMieOpticalDepth  = g_MediaParams.f4MieExtinctionCoeff.rgb      * f2NetParticleDensityFromCam.y;\n"
"    f3Extinction = exp( -(f3TotalRlghOpticalDepth + f3TotalMieOpticalDepth) );\n"
"\n"
"    // Apply phase function\n"
"    // Note that cosTheta = dot(DirOnCamera, LightDir) = dot(ViewDir, DirOnLight) because\n"
"    // DirOnCamera = -ViewDir and LightDir = -DirOnLight\n"
"    float cosTheta = dot(f3ViewDir, f3DirOnLight);\n"
"    ApplyPhaseFunctions(f3RayleighInscattering, f3MieInscattering, cosTheta);\n"
"\n"
"    f3Inscattering = f3RayleighInscattering + f3MieInscattering;\n"
"}\n"
"\n"
"\n"
"// BACKUP\n"
"//\n"
"// To better distribute samples for numerical integration I tried used the following simplifying assumptions:\n"
"// - Uniform media density\n"
"// - No light extinction\n"
"// Under these assumptions the total amount of light inscattered from 0 to distance t will be computed by\n"
"//\n"
"//      I(t) = Integral(exp(-a*x), 0 -> t) = 1/a * (1 - exp(-a*t))\n"
"//\n"
"// We want to find sampling that produces even distribution for function I(t):\n"
"//\n"
"//      I(tn) = I(D) * n / N, where D is the total distance and N is the total sample count.\n"
"//\n"
"// From here:\n"
"//\n"
"//      1/a * (1 - exp(-a*tn)) = I(D) * n / N\n"
"//      1 - exp(-a*tn) = a * I(D) * n / N\n"
"//      exp(-a*tn) = 1 - a * I(D) * n / N\n"
"//      tn = -1/a * log(1 - a * I(D) * n / N) = -1/a * log(1 - n/N * (1 - exp(-a*D)))\n"
"//\n"
"// This unfortunately did not work because homogeneous media assumption is inappropriate:\n"
"// -  when a = 1e-5, the distribution is basically linear\n"
"// -  when a = 1e-4, the first N-1 samples cover first 50% distance, while the last one covers the remaining 50% distance\n"
